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ABSTRACT 

We investigate the structure of the dark matter halo formed in the cold dark 
matter scenario using A^-body simulations. We simulated 12 halos with the 
mass of 6.6 x IO^^Mq to 8.0 x lO^^M©. In almost all runs, the halos have density 
cusps proportional to r~^'^ developed at the center, which is consistent with the 
results of recent high-resolution calculations. The density structure evolves in a 
self-similar way, and is universal in the sense that it is independent of the halo 
mass and initial random realization of density fluctuation. The density profile is 
in good agreement with the profile proposed by Moore et al. (1999), which has 
central slope proportional to r~^'^ and outer slope proportional to r~^. The halo 
grows through repeated accretion of diffuse smaller halos. We argue that the 
cusp is understood as a convergence slope for the accretion of tidally disrupted 
matter. 



Subject headings: cosmology: theory — dark matter — galaxies:kinematics and 
dynamics — galaxies:formation — method: numerical 



1. Introduction 

In standard cosmological pictures, such as the cold dark matter cosmology, dark matter 
halos are considered to be formed in a hierarchical way; smaller halos first formed from 
initial density fluctuations and they merged with each other to become larger halos. In 
reality, the formation process of dark matter halo is rather complicated, since a variety of 
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processes, such as merging between halos of various sizes and tidal disruption of small halos 
(satellite) proceed simultaneously. 

One of the most influential works on the dark matter halo is the "finding" of the 
universal profile by Navarro, Prenk, and White (1996, 1997, hereafter NFW), though there 
were many analytical and numerical studies before NFW (see NFW (1996) or Bertschinger 
(1998) for reviews). NFW performed A^-body simulations of the halo formation and found 
that the profile of dark matter halo can be fitted by a simple formula 



(r/rs)(l + r/rs)^' 

where po is a characteristic density and Vg is a scale radius. They also argued that the 
profile has the same shape, independent of the halo mass, the initial density fluctuation 
spectrum or the value of the cosmological parameters. It should be noted that, before 
NFW, Dubinski and Carlberg (1991) also found in their high-resolution simulation the halo 
can be well fitted by Hernquist (1990) profile . 

Many studies on the NFW "universal profile", both numerical and analytical, were 
done after their proposal. Many A^-body simulations whose resolution are similar to those 
of NFW were performed and results similar to NFW were obtained (Cole and Lacy 1996, 
Tormen, Bouchet and White 1996, Brainerd, Goldberg, and Villumsen 1998, Thomas et al. 
1998, Okamoto and Habe 1999, Huss, Jain and Steinmetz 1999, Kravtov et al. 1998, Jing 
2000). Analytical and semi- analytical studies to explain the NFW universal profile were 
also done (Evans and Collet 1997, Syer and White 1998, Avila-Reese, Firmani, Hernandez 

1998, Nusser and Sheth 1999, Kull 1999, Heriksen and Widrow 1999, Yano and Gouda 

1999, Bullock et al. 1999, Subramanian, Cen, Ostriker 2000, Lokas 2000). The clear 
understanding for the NFW profile, however, has not yet been given. One of the reasons 
why a clear understanding has not been established might be that all these studies were 
trying to answer a wrong question. 

Our previous study (Fukushige and Makino 1997, hereafter FM97) showed that density 
profile obtained by high-resolution iV-body simulation is different from the NFW universal 
profile. We performed simulations with 768k particles, while previous studies employed 
~ 20k. We found that the galaxy-sized halo has a cusp steeper than p oc r~^. 

This disagreement with the NFW universal profile was confirmed by other high- 
resolution simulations. Moore et al. (1998, 1999) and Ghigna et al. (2000) performed 
simulations with up to 4M particles and obtained the results similar to ours. They found 
that cluster-sized halos also have cusps steeper than the NFW profile and they proposed 
the modified universal profile, p — po/[('^/'^s)^'^(l + ('^/'^s)''^'^)]- On the other hand, Jing and 
Suto (2000) found that the density profile of dark matter is not universal. They performed 
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a series of A'"-body simulations and concluded that the power of the cusp depends on mass. 
It varies from —1.5 for galaxy mass halo to —1.1 for cluster mass halo. 

In this paper, we again investigate the structure of dark matter halos using A^-body 
simulation. We performed A^-body simulations of formation of 12 dark matter halos with 
masses 6.6 x IO^^Mq to 8.0 x IO^^Mq, using a special- purpose computer GRAPE-5 (Kawai 
et al. 2000) and Barnes-Hut treecode. In section 2, we describe the model of our A'^-body 
simulation. In section 3, we present the results of simulation. Section 4 is for summary and 
section 5 is for discussion. 



2. Simulations Models 

We performed in total 12 runs on 4 different mass scales, which are summarized 
in Table 1. Initial conditions were constructed in a way similar to that in FM97. We 
assigned initial positions and velocities to particles in a spherical region with a radius of R 
Mpc surrounding a density peak selected from a unconstrained discrete realization of the 
standard CDM model = 50km/s/Mpc, Q = 1 and (jg = 0.7). The peak was chosen from 
an i?box Mpc cube using a density field smoothed by a Gaussian filter of radius (i?box/2) 
Mpc. The values of R and i?box in the comoving flame are summarized in Table 1. In order 
to generate the discrete reahzation of the CDM model we used the COSMICS package. 



Table 1: 


Simulation 


Models 












Run 


R (Mpc) 


i^box (Mpc) 


m (Mo) 


e (kpc) 


At (yr) 


^-start 


-^end 


16M0 


12.8 


32 


3.0 X 10^ 


0.56 


1.6 X 10^ 


18.8 


0.0 


16M1 


16 


32 


6.0 X 10^ 


0.56 


1.5 X 10^ 


18.5 


0.0 


16M2 


16 


32 


6.1 X 10^ 


0.56 


1.6 X 10^ 


20.4 


0.0 


8M0 


6.4 


16 


3.7 X 10^ 


0.28 


7.9 X 10^ 


22.3 


0.58 


8M1 


8 


16 


7.6 X 10^ 


0.28 


7.5 X 10^ 


22.2 


0.63 


8M2 


8 


16 


7.6 X 10' 


0.28 


7.8 X 10^ 


23.9 


0.59 


4M0 


3.2 


8 


4.7 X 10^ 


0.14 


7.8 X 10^ 


25.9 


1.6 


4M1 


4 


8 


9.5 X 10^ 


0.14 


7.4 X 10^ 


25.9 


1.6 


4M2 


4 


8 


9.5 X 10^ 


0.14 


7.6 X 10^ 


27.4 


1.2 


2M0 


1.6 


4 


5.9 X 10^ 


0.07 


3.8 X 10^ 


29.7 


2.1 


2M1 


2 


4 


1.2 X 10*^ 


0.07 


3.6 X 10^ 


29.7 


2.2 


2M2 


2 


4 


1.2 X 10^ 


0.07 


3.8 X 10^ 


30.9 


1.8 



We followed evolution of the density peak by N-hody simulation. We added the local 
Hubble flow and integrated the orbits directly in the physical space. We used the Plummer 
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softened potential with the softening length constant in physical space, and used a leap-flog 
integrator with shared and constant timestep. In Table 1 we summarized the individual 
particle mass, m, softening length, e, timestep size. At, and starting and ending redshift, 
2;start and Zend- The masses of particles are equal and the total number of particles for each 
simulation is (2.0 - 2.1) x 10^. 

We determined the radius R Mpc for Run 16M{0,1,2} using trial runs with smaller 
number of particles, so that all particles lying inside of r2oo at z^nd are included. Here, the 
radius r2oo is defined as the radius of the sphere in which the mean density p is equal to 
200pcrit, where pcrit is the critical density. We did not include tidal effects from outside the 
R Mpc sphere. The region of Runs 8M{0,1,2}, 4M{0,1,2}, and 2M{0,1,2} is 1/2, 1/4, and 
1/8 of the size of 16M{0,1,2} and mass resolution are x2, x4, and x8, respectively. The 
ending redshift Zend for Runs 8M{0,1,2}, 4M{0,1,2}, and 2M{0,1,2} is determined so that 
the truncation outside the sphere did not infiuence the profile around r2oo- 

The number before M in a run name (ex. 8 for Run 8M1) indicates the length of the 
simulation box, i?box/2, in Mpc. The number after M in the run name (ex. 1 for Run 8M1) 
identifies the index for the random number seed used to generate initial density field. For 
example. Runs {16,8,4,2}M0 are series of runs in which the phases of initial density waves 
are same and the amplitude of the wave are different. 

For the force calculation, we used the Barnes-Hut tree code {6 = 0.75, Barnes and 
Hut 1986, Barnes 1990, Makino 1991) implemented on GRAPE-5 (Kawai et al. 2000), a 
special-purpose computer designed to accelerate A^-body simulations. Using the tree code 
on two GRAPE-5 boards and a workstation whose CPU is 21264/677MHz Alpha chip, one 
timestep took 21 seconds. The total number of timesteps is about 6000 — 8000. Therefore, 
we can complete one run in 35-50 CPU hours. 

3. Result 
3.1. Snapshot 

Figure shows the particle distributions for Run 16M0 at 16 different redshifts. For 
these plots, we shifted the origin of coordinates to the position of the potential minimum so 
that the largest halo is at the center of the panel. Figure ^ shows the particle distribution 
for Runs {16,8,4,2}M1 together. The phases of waves for initial density field are the same 
for all these runs and only the amplitude are different. In Table 2, we summarized the 
radius r2oo, the mass M200 and number of particles A'200 within r2oo, at Zend- 
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3.2. Accuracy Criteria 

In this study, we plot the density only for the radii unaffected by numerical artifacts. 
We used the following two criteria: (1) tj.ei{r)/t > 3 and (2) tdy{r)/At > 40. We obtained 
the criteria (1) and (2) experimentally and the details of the experiments are discussed in 
sections |3.2.1| and |3.2.2| . We plot the densities only if both criteria are satisfied. Here, trei{r) 



is the local two-body relaxation time defined by, 

0.065f3 



''rel 



G^pm ln(l/e) ' 

(cf. Spitzer 1987) and tdy('^) is the local dynamical time defined by 



(2) 



= (Gp)-'/' (3) 

where p is the average density within radius r. Using these criteria we judge whether the 
density profile is unaffected by numerical artifacts due to the two-body relaxation (1) and 
the step size for the time integration (2). 

The inner limit of radii where the density is correctly calculated in our simulations is 
(0.01 — 0.02)r2oo at the final redshifts. In most cases, the criterion (1) for the two-body 
relaxation determines the limit radius for reliability. 



3.2.1. Criterion for two-body relaxation 

In this subsection, we evaluate the criterion, trci(^)/i > 3, to distinguish the numerical 
artifact due to the two-body relaxation. In order to see the effect of two-body relaxation, 
we calculated the same model as Run 16M0 but with several different values for total 
number of particles (A^) and the softening size (e). In Figure ^, we plot the final average 
density profiles p for three simulations with N/4, N/16, and A^/16 and Ae, where N and e 
mean values used in Run 16M0. Otherwise stated, we used the same simulation parameter 
as in Run 16M0. We can see that the central density depends on rather strongly, and is 
lower for lower number of particles. For the same value of N, larger softening has small but 
clear effect of increasing the central density. 

Figure ^ shows the ratio, p/pref, where pref is the averaged density of the reference run 
in which the effect of the two-body relaxation is smallest. Here, we used Run 16M0 as 
the reference run. In Figure ^, we show the ratio p/ prcf plotted as a function of the ratio 
of the local two-body relaxation time t^ei{f) defined by (0) to simulation period t. Note 
that trei('") is monotonous and increasing function of r. Therefore, smaller tj-^ilr) means 
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smaller r. We can see that the density tends to go below the reference value if t < ?,t^e[{r). 
The difference between the reference run and runs with smaller number of particles is 
insignificant if t < ?)t^e\{f). From this result, we adapt tj.ci{r)/t > 3 as the criterion for the 
two-body relaxation. 

3.2.2. Criterion for time integration 

In this subsection, we evaluate the criterion, tdy{r)/At > 40, to distinguish the 
numerical artifact due to large step size of time integration. In order to see the effect of 
large step size, we calculated the same model as Run 16M0 but with larger time step size 
(At). In Figure |^, we plot the final average density profiles p for three simulations with 
4At, 8At, and 16At, where At means values used in Run 16M0. Otherwise stated, we 
used the same simulation parameter as in Run 16M0. We can see that the central density 
depends on At rather strongly and is lower for larger At. Figure |^ shows the ratio, p/prei- 
Here, we used Run 16M0 as the reference run. 

In Figure ^, we show the ratio p/pref plotted as a function of ratio of the local 
dynamical time t^y defined by (H) to the time step size At. Note that tdyl'^) is monotonous 
and increasing function of r. Therefore, smaller t^j{r) means smaller r. We can see that the 
density tends to go below the reference value if t^yir) < 40At. From this result, we adapt 
tdy{f')/At > 40 as the criterion for time integration. 

Note that the number 40 is applicable only to the integration scheme we used: the 
leapfrog scheme integrated in physical space with constant time step size. This scheme has 
good characteristics such as the time-reversibility and symplecticity. The number 40 should 
increase when variable stepsize is used or the system is integrated in comoving space (where 
acceleration depends on velocity). 

As a result of adapting this criterion, the total number of time steps to integrate in 
Hubble time becomes about 8000 in our simulation. The number is a little smaller than 
that reported in previous simulation (ex. ~ 50000, Moore et al. 1999). We also calculated 
the same model as Run 16M0 but with 4 times many time steps. We confirmed that the 
density profile, shape, and anisotropy parameter do not change, and the number of timestep 
8000 is enough. In figure ^ we show the density profile. 
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3.2.3. Other numerical effects 

The potential softening also affects the density profile. In order to see the effect of the 
potential softening, we calculated the same model as Run 16M0 but with different softening 
length [e). In Figure ^ we plot the final average density profiles p for six models with e = 
0.18, 0.56, 1.7, 5, 15, 30 kpc. Except for the softening length, we used the same simulation 
parameter as in Run 16M0. In Figure ^ we can see that the central density is lower for 
both smaller and larger e. The former is because the two-body relaxation effect is stronger 
and the time integration is less accurate for smaller e. The latter is because the potential 
softening itself affects the density structure for larger e. The potential softening, therefore, 
should be set in the intermediate range. 

The potential softening for Run 16M0 {e =0.56kpc) is in the intermediate range, 
though it is not optimal. It does not affect the density profile outside of 20 kpc, which is 
the critical radius defined by the accuracy criterion (1) and (2). The ratio of the softening 
length to the critical radius in other runs are similar to that in Run 16M0. 

We made sure that the accuracy of the BH tree code did not influence the structure in 
the range where the above criteria are satisfied, by re-simulating the same initial model as 
used in FM97, in which the direct summation is used. We found no systematic difference in 
the results. Therefore, the accuracy of BH tree-code is okay. 



3.3. Density Profiles 

Figure |I^ shows the evolution of density profiles for Run 16M0. The position of the 
center of the halo was determined using the potential minimum and the density is averaged 



over each spherical shell whose width is log]^Q(Ar) = 0.0125. Figure 11 is for Run 2M0. 
For the illustrative purpose, the densities are shifted vertically. Figure |12] show the density 
profiles at Zend for all Runs. 

In all runs we can see the central density cusps approximately proportional to r^^'^. 
In other words, the power of the cusp is —1.5 and is independent of halo mass, which is 
consistent with the result of Moore et al.(1999). In the outer region, the density profiles are 
very similar for all runs. The dependence of the power-law index of the inner cusp on the 
halo mass observed by Jing and Suto (2000) was not reproduced in our simulations. Even 
if we take into account the run-to-run variation, the dependence on mass in our results 
is in the opposite direction compared to that of Jing and Suto (2000). In the following 
subsections, we discuss the self-similar growth of the halo (section p.4|) , the universahty of 
the profile (section p.5| ), and the mechanism for the self-similar growth (section |3.6| ). 
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3.4. Self-Similar Evolution 



Figure shows the growth of the halo, without the vertical shift. In this figure it is 
clear that the halo grows in a self-similar way, keeping the density of the central cusp region 
constant. 

If the evolution is self-similar, we can write the density as 

p(r,t) = Pt(^)P*(^*) (4) 
r, = r/rt(M) (5) 

Here, we write Pt(M) and r^{M) as a function of the mass of the halo M, instead of the 
time. The self-similar profile itself should have the central cusp of p*(r*) oc r". The actual 
profile at the cusp region satisfies p(r) = Cr", with C constant in time. Therefore, p| and 
Tj. should satisfy oc r". If we write p-^ and as function of M, we have 

Pm = Poo(^) (6) 



00 



where roo, PoO) ^oo are constants and n is the power-law index of the cusp given by p oc r". 



This self- similarity is illustrated in Figure 0. If we set n = —1.5 from the simulations, we 
obtain 

Pt(M) oc M-^ (8) 
rt(M) oc Mi. (9) 



In Figure |T5|, we plot defined through equations (H)-® as a function of Here, 
we took Moo = roo = 0.2 Mpc, and poo = 7 x IO'^Mq/pc^. We plot four density 

profiles at different values of the redshift z. We set n = —1.5 for all runs and use M200 as the 
total mass. In this figure, we can see that the density profiles of the same halo at different 
times show very good agreement to each other, which means the density structure evolves 
self-similarly, though in outer region a degree of overlapping becomes worse. Therefore, 
Figure |T3| demonstrates that our assumption of the self- similarity is justified. 



3.5. Universality 

In this subsection we discuss the universality of the density profile. Using a 
non-dimensional free parameter 6, we define new non-dimensional variables expressed as 

p** = p*6^^ (10) 
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r** = r*(5» (11) 

Figure |l6| shows p** r*^, of all Runs at z = Zcnd- The values of 6^^ are 1.0, 0.4, and 0.6 for 
Run 16M{0,1,2}, 2.5, 1.0 and 3.0 for Run 8M{0,1,2}, 10.0, 3.0 and 6.0 for Run 4M{0,1,2}, 
and 35.0, 12.0 and 30.0 for Run 2M{0,1,2}. We can see that the 12 density structures agree 
very nicely, which means they are universal. In principle, any scaling on the r-p plane can 
be expressed using two parameters. We used the total mass M as one of two parameters 



to express the self-similarity discussed in section |3.4| . The parameter S corresponds to the 
other freedom. The value of S is considered to reflect an amplitude of the density fluctuation 
at the collapse. 

We attempted to fit the density structure to several profiles proposed in earlier studies. 
Figure ^ and |T8| show the profile proposed by Moore et al. (1999) and by NFW. The 
function forms are given by p** = r~^'^{l + rlf)~^ and by p** = 10(0.5r^,^,)~^(l + 0.5r**)~^, 
respectively. Our simulation results agree with the profile proposed by Moore et al.(1999) 
very well, while the agreement with the NFW profile is not so good. 

In summary, the density structure of simulated halos is well expressed by 

77^^ ^ (r/ro)i-5[l + (r/ro)i-5] ^^^fTT^]) ^^^^ 

where 



Po = 7x10-^-5 -—— (Mo/pc^) (13) 




^0 = 0-2 ■ ^"M Y^n]^^ I (Mpc) (14) 

where again M is the total mass of halos, 5 is a non-dimensional free parameter. The free 
parameter d is constant during evolution of a halo. 



3.6. Formation Process of The Central Cusp 

In this subsection, we show that the central cusp grows through the accretion of the 
disrupted smaller halos by a larger halo. In the bottom-up structure formation a typical 
halo grows through repeated merging of smaller halos. In the CDM hierarchical clustering, 
the larger halo typically has a denser central region than the smaller halo has. Therefore, 
the smaller halo is disrupted by the tidal field of the larger halo and the matter from the 
smaller halo is scattered around, when two halos merge. On the other hand, the central 
region of the larger halo survives the merging process more or less intact. 
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in Figure 13, we can clearly see that the cusp grows outward without changing the 
inner part. In Figure |19| we show density profiles of halos which will merge to the largest 
halo, for Runs 16M0 and 2M0. We plot 6 halos with more than 1000 particles which are 
nearest from the potential minimum, together with the largest halo. It is clear that the 
central halo has the highest density. 

The reason why larger halos have higher density can be understood as follows. Let us 
consider the peaks of the fluctuation whose characteristic scale is A. If the total density field 
is only composed of the fluctuation whose scale is A, the peaks will collapse to halos with 
similar density almost simultaneously. Actually, there are contributions from fluctuations 
whose scale is larger than A, too. A peak in high-density background would collapse to 
a halo with higher density than peaks in low-density background, simply because of the 
difference in the background density. Later, the "background", which is just a density peak 
of longer wavelength, would collapse. During this collapse, however, the high-density peak 
collapsed earlier is not affected. Therefore, larger halo tend to have higher central density 
than smaller halos. 

In Figure ^ we show one- dimensional trajectory of the particles for Run 16M0. We 
randomly select 10 particles from the particles whose distances from the center of the halo 
at the end of the simulation are 0.02-0.03, 0.1-0.2, and 1-2 Mpc Figure ^ shows that a large 
fraction of particles in the inner region settles there early, and those in the outer region 
tend to fall later. In other words. Figure ^ shows that the formation process discussed in 
the above actually takes place. 

The cusp with the slope of —1.5 seems to be a "fixed point" or a "convergence point" 
for the growth of the halo through accretion of diffuse and small halos. Once the cusp with 
the slope of —1.5 forms, the density in the r~^'^ cusp remains unchanged and the disrupted 
matter is accreted outside the r~^'^ cusp, which is clearly shown in Figure p!3[ 

Moreover, the power index of —1.5 seems to be a universal feature independent of the 
form of initial power spectrum. The high-resolution simulations presented in this paper and 
those by Moore et al. (1999) show that the power of the cusp is —1.5, independent of the 
mass scale. A preliminary result of our another simulation from the initial power spectrum 
of P{k) oc fc"^ '', which is shallower than that at cluster scale for standard CDM model, also 
show that the power of the cusp is around —1.5. 

Currently, we do not have a clear explanation why the slope of the cusp is —1.5 when 
it forms through the accretion of disrupted small halos. We will discuss this topic more 
comprehensively elsewhere. 
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3.7. Origin of The Outer Profile 



Figure ^ shows the distribution of particles on the r-f ^ plane, where r is distance from 
the center and is the radial velocity, at 16 different redshifts for Run 16M0. We can see 
that the outer region consists of two component. The first component is infalling matters 
which is visible as thick stream of particles in the right-lower region of each panels. The 
vertical spreads visible in this stream are infalling smaller halos. The second component 
is the more scattered particles with nearly zero average velocity. As one can see from the 
time evolution, these stars gained energy in the central region when small halos accreted on 



the central halo. In Figure 22 we show density profiles of scattered and infalling particles 



separately. We separated two components by defining appropriate boundary in Figure 21 



At around r2005 contribution of two components to the total profile are of the same order. 

The profile in the outer region exhibits large fiuctuations. The merging events occur 
intermittently, and the amount of scattered matter depends on earlier merging events. 
Nevertheless, the density profile fits the profile which is asymptotically proportional to r~^, 
as shown in the previous section. 

Consequently, in the outer region orbits of particles shows strong radial anisotropy. 



Figure |23| shows anisotropy in velocity distribution of the profile together with simulation 
results for Run 16M{0,1,2} and 2M{0,1,2}. The anisotropy are expressed by the anisotropy 
parameter /3, defined as 

^-l-f. (15) 

where and (f^) are mean tangential and radial velocity dispersion. In this definition, 
(3 = means the velocity distribution is isotropic and (3 = 1 means completely radial. 



4. Conclusion 

We performed A^-body simulations of dark matter halo formation in the standard CDM 
model. We simulates 12 halos whose mass range is 6.6 x lO^^M© to 8.0 x 10^^ Mq. We 
introduced the accuracy criteria to guarantee that numerical artifact due to the two-body 
relaxation and the time integration do not affect the result, and obtained the density profile 
which is free from numerical artifact down to the radii (0.01 — 0.02)r2oo- 

Our main conclusions are: 
(1) In all runs, the final halos have density cusps proportional to r~^'^. 
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(2) The density profile evolves self-similar ly. 

(3) The density profile is universal, independent of the halo mass, initial random 
reahzation of density fluctuation and the redshift. The density structure is in good 
agreement with the profile proposed by Moore et al. (1999). 

(4) We found that the central cusp grows through the disruption and accretion process of 
diffuse smaller halos. The slope of the central cusp seems to be a fixed point for the 
growth of the halo through accretion of tidally disrupted matter. 

5. Discussion 

Here we discuss the relation between our results and those of the previous studies. 

We obtained steeper inner cusp than that obtained by NFW, which was already 
found in high-resolution simulations (FM97, Moore et al 1999, Ghigna et al.2000, Jing 
and Suto 2000). The reason for this disagreement is that in low- resolution simulations the 
central cusp is smoothed out by the two-body relaxation. If the cusp is shallower than —2, 
the velocity dispersion decreases inward. The energy flows inward due to the two-body 
relaxation and the central region expands, which is called the gravothermal expansion 
(Hachisu et al. 1978, Quinlan 1996, Heggie, Inagaki, McMillan 1994, Endo, Fukushige, 
Makino 1997). Therefore, the density in the cusp decreases and the cusp becomes shallower. 
Using the relations irei ~ '^^/{p^)-i P ~ r~^-^, and v ~ r^^'^^ and our simulation results, we 
can estimate the lower limit of the radius where the structure is free from the two-body 
relaxation effect as ~ 0.01r2oo(A^/10^)"°-^^(po/2.7 x lQ-^MQ/^c^f-^^{ao/l2>m'kmls)-^-^^ for 
a cluster-sized halo at the present epoch, where ctq is velocity dispersion at the scale radius 
tq. For simulations with = 10^ and 10^ within r2005 the limits are estimated as ~ 0.08r2oo 
and ~ 0.03r200; respectively. Therefore, simulations with ~ 20k particles the central cusp 
would be become signiflcantly shallower due to relaxation. 

We could not reproduce the dependence of the slope observed by Jing and Suto 
(2000). The difference could be also due to the smoothing by two-body relaxation in 
their cluster-sized halos. In this paper, we show that the density proflic within ~ 0.01r2oo 
smoothed by the two-body relaxation. The density at 0.01r2oo and the mass resolution in 
their cluster-sized halo are similar to those in ours. The density proflle in their simulations 
within 0.01r2oo, at which the proflle begins to depart from r~^'^ inward, could be affected 
by the two-body relaxation. 

Our result is in good agreement with results of simulations by Moore et al. (1999), in 
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which the tidal field was included. This agreement suggests that the neglect of the tidal 
field in our present study and in FM97 hardly affects the density profile. The formation 
mechanism we discussed in this paper also suggests that the tidal field due to the mass 
outside the simulation sphere is not crucial. 

Moore et al. (1999) argued that the merging process is not related with the structure. 
They simulated the halo formation using a power spectrum with a cutoff to suppress 
merging event in smaller scale and showed the profile does not change. However, in this 
simulation, several merging event took place since the cutoff wave length is rather short. 
Therefore, their conclusion that merging is not important is not really supported by their 
simulation. According to our explanation, several merging events where the central large 
halo swallows smaller halos determines the structure of the halo. Such events did took place 
in Moore et al.'s simulations. 

As discussed in the Introduction, there are a lot of analytical and semi-analytical 
studies to explain the density structure. However, no study succeeded to explain the 
universality satisfactory. This is because none of them is based on the formation and growth 
process of the halo as discussed in this paper. Syer and White (1998) discussed that the 
cusp is a convergence point of disrupting and sinking of satellite, and that the convergence 
slope depends on the initial power spectrum. In their model, some of smaller halos were 
assumed to sink down to the center of the larger halo, when two halos merged. However, 
as we discussed, smaller halos are always disrupted and never sink down to the center, 
because the smaller halo is always less dense. Evans and Collet (1997) argued the cusp 
is a steady-state solution of Fokker-Plank equation. The solution is derived by assuming 
that many small clumps within a large halo evolve by the two-body relaxation. In our 
simulations, there are no such small clumps, since they are disrupted before they reaches 
the center. 

We are grateful to Atsushi Kawai, for his help in preparing the hardware and software 
environment of the GRAPE-5 system, and to Yasushi Suto and Yoko Funato for many 
helpful discussions. To generate initial condition, we used the COSMICS package developed 
by Edmund Bertschinger, to whom we express our thanks. A part of numerical computations 
was carried out on the GRAPE system at ADAC (the Astronomical Data Analysis Center) 
of the National Astronomical Observatory, Japan. This research was partially supported by 
the Research for the Future Program of Japan Society for the Promotion of Science, grant 
no. JSPS-RFTP 97P01102. 
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Fig. 1. — Snapshots from Run 16M0 at 16 different redshifts indicated by the number at 
lower-left of each panel. The length of the side for each panel is equal to 4r2oo at Zend- 
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Fig. 2. — Snapshots from Runs {16,8,4,2}M1 at the final redshifts. The length of the side 
for each panel is equal to 4r2oo ^end- 
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Fig. 3. — Averaged density profiles for Run 16M0 and tlie same model as Run 16M0 but witli 
and iV/16 and 4e, wfiere N and e are tfie number of particles and softening 
parameter for Run 16M0. The model with Ae is plotted in the dashed line. The unit of 
density is Mq/pc^. 




.01 .03 .1 .3 1 

r(Mpc) 



Fig. 4. — The ratio of the averaged density for the same models as Fig ^ to that for the 
reference run (Run 16M0). 



- 19 - 




Fig. 6. — Averaged density profile for Run 16M0 in tlie tliick line and the same model as 
Run 16M0 but with 4, 8, and 16 times larger and 4 times smaller time stepsize in the thin 
line. The unit of density is Mq/pc^. 



- 20 - 




.01 .03 .1 .3 1 

r(Mpc) 



Fig. 7. — The ratio of the averaged density for the same model as Fig ^ for the reference run 
(Run 16M0). 
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Fig. 8. — Same as Figure ^ but as a function of the ratio of the local dynamical time to 
time step size, t^y{r)/ At. The arrow indicates the accuracy criterion 40. 
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Fig. 9. — Averaged density profiles for Run 16M0 {e = 0.56kpc) and tlie same model as Run 
16M0 with e =0.18, 1.7, 5, 15, 30 kpc, where e is the softening length. The unit of density 
is M0/pc^. The numbers beside the lines indicates the softening length. The profiles for 
Run 16M0 and the models with e — 1.7kpc are plotted in the thick lines. The dashed hues 
indicates the density profile proportional to r~^"^. The arrow indicates the critical radius 
defined by the accuracy criteria (1) and (2). 
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Fig. 10. — Evolution of density profile of the halo for Run 16M0. The unit of density is 
Mq/pc^. The profiles are vertically shifted downward by 7, 6, 1, dex from the bottom to 
the top. The arrows indicates r2oo- The numbers near the dashed lines indicate the power 
index of those lines. The numbers on the left of the profiles indicate the redshift. 
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Fig. 12. — Density profile of the fialo for all Runs at Zf^-^A- The unit of density is Mg/pc^. 
The profile for Run 16M{0,1,2}, 8M{0,1,2}, and 4M{0,1,2} are vertically shifted upward 
by 1.5, 1, 0.5 dex, respectively. The arrows indicate r2oo- The thin sohd and dashed hues 
indicate the densities proportional to T~^^^ and . 
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Table 2: Halo properties at z = Zend 



Run 


M200 (Mq) 


r2oo (Mpc) 


^200 


16M0 


2.6 X 


1.7 


873170 


16M1 


7.8 X 10^4 


2.4 


1279383 


16M2 


8.0 X 10^^ 


2.4 


1322351 


8M0 


2.8 X 10^3 


0.48 


745735 


8M1 


9.0 X 10^3 


0.72 


1186162 


8M2 


7.7 X 10^^ 


0.70 


1015454 


4M0 


2.7 X 10^2 


0.13 


559563 


4M1 


8.0 X 10^2 


0.20 


846301 


4M2 


6.6 X 10^2 


0.22 


697504 


2M0 


6.6 X 10^^ 


0.062 


643151 


2M1 


1.1 X 10^^ 


0.085 


957365 


2M2 


1.0 X 10^2 


0.096 


923545 
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Fig. 13. — Evolution of the density profile of the halo for Runs 16M0 and 2M0 with no 
vertical shift. The unit of density is Mq/pc^. 



Fig. 14. — Self-Similar Evolution. 



-28- 




Fig. 15. — Self-similar evolution of the density profile. The scaled densities are plotted as 
a function of the scaled radius r*. The profiles for Run 8M{0,1,2}, 4M{0,1,2} and 2M{0,1,2} 
are vertically shifted downward by 0.5, 1, 1.5 dex and horizontally shifted to the left by 0.5, 
1, 1.5 dex, respectively. 
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Fig. 17. — Fitting of the density structure by the profile proposed by Moore et al. (1999) 
(solid curve). The dashed line indicates p oc r~^^^. 
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Fig. 18. — Fitting of the density structure by the profile proposed by Navarro, Frenk, White 
(1996,1997). 
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Fig. 19. — Density profiles of smaller halos which are going to merge to the larger halo. The 
thick curve indicates the larger halo. 




Fig. 20. — One-dimensional trajectories of randomly selected 10 particles for Run 16M0. 
The distance from the center of halo is plotted as a function of time, together with the 
density profile at Zend- The region indicated by thick curve in density profile indicates where 
the selected particles exists at 2;end- 
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Fig. 21. — Distribution of particle in r — Vr plane at 16 different redshifts for A'"-body 
simulation Run 16M0. The number at lower-left corner of each panel indicates the redshift. 
The width of each panel is equal to 2r2oo at 2;end- The range of radial velocity is from —2031 
km/s to 2031 km/s. 
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Fig. 22. — Density profile of the halo for Run 16M0 as a superposition of two components. 
The dashed curves indicate that for the infaUing matter and the thick curves indicate the 
scattered matter. The thin curve is total profile. 
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Fig. 23. — Anisotropy parameter, (3 = 1 — v^/v^ as a. function of radius at 2; = ^end for Runs 
16M{0,l,2}(solid curve) and 2M{0,1,2} (dashed curves). The star, square, and triangle 
symbol indicate P for {16,2}M0, {16,2}M1, and {16,2}M2, respectively. 



